
####################################################
####          Replication Code for              ####
####     Bayesian versus maximum likelihood     ####
####      estimation of treatment effects       ####
####          in bivariate probit               ####  
####       instrumentalvariable models          ####    
####################################################
####                                            ####  
####                4/ 18/2017                  ####         
####     this file loads the STATA ML results   ####
####      and calculates combines estimates     ####
####################################################


###load results
biprobit <- read.xls("results_biprobit.xlsx", header = F)
names(biprobit) = c("perry", "ihad", "star")



sumVar = ((0.8275299^2) + (1.031282^2)+(.9214786^2))

sumVar = ((biprobit$perry[8]^2) + (biprobit$ihad[8]^2)+(biprobit$star[8]^2))

wPerry = (biprobit$perry[8]^2 )/sumVar
Perry_coef =  biprobit$perry[7] 

wIhad = (biprobit$ihad[8]^2)/sumVar
ihad_coef = biprobit$ihad[7]

wStar = (biprobit$star[8]^2)/sumVar
star_coef = biprobit$star[7]

#weights sum to one:
sum(c(wStar, wIhad, wPerry))

### combined coef:
combined_coef = ((star_coef*wStar) + (ihad_coef*wIhad)  +  (Perry_coef*wPerry))

## combined se:
combined_SE = sqrt(sumVar)/3

### here is what was done in the paper:
#wStar2<-(0.880)
#wIhad2<-(1.070)
#wPerry2<-(0.888)
#this<-((0.880^2) + (1.070^2)+(0.888^2))
#((1.948*wStar2) + (1.051*wIhad2)  +  (1.084*wPerry2))/this
#Gives 1.4038 
## i.e. they did not square the standard error
##((1.948*wStar2^2) + (1.051*wIhad2^2)  +  (1.084*wPerry2^2))/this
## This is what they should have gotten





